Pressure-induced diamond to (3-tin transition in bulk silicon: 
a near-exact quantum Monte Carlo study 
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The pressure-induced structural phase transition from diamond to /3-tin in silicon is an excellent test for 
theoretical total energy methods. The transition pressure provides a sensitive measure of small relative energy 
changes between the two phases (one a semiconductor and the other a semimetal). Experimentally, the transition 
pressure is well characterized. Density-functional results have been unsatisfactory. Even the generally much 
more accurate diffusion Monte Carlo method has shown a noticeable fixed-node error. We use the recently 
developed phaseless auxiliary-field quantum Monte Carlo (AFQMC) method to calculate the relative energy 
differences in the two phases. In this method, all but the error due to the phaseless constraint can be controlled 
systematically and driven to zero. In both structural phases we were able to benchmark the error of the phaseless 
constraint by carrying out exact unconstrained AFQMC calculations for small supercells. Comparison between 
the two shows that the systematic error in the absolute total energies due to the phaseless constraint is well within 
0.5m£'h/atom. Consistent with these internal benchmarks, the transition pressure obtained by the phaseless 
AFQMC from large supercells is in very good agreement with experiment. 

PACS numbers: 64.70.K-, 71.15.-m, 61.50.Ks, 71.15.Nc. 



I. INTRODUCTION 

Theoretical and computational treatment of the effects of 
electron correlations remains a significant challenge. De- 
spite decades of effort invested into solving the Schroedinger 
equation (by independent-particle, mean-field and perturba- 
tive methods), there are still major difficulties in predicting 
and explaining many phenomena related to bonding, cohe- 
sion, optical properties, magnetic orderings, superconductiv- 
ity and other quantum effects. The pressure-induced structural 
phase transition in silicon from diamond to /3-tmi is an excel- 
lent test for theoretical total energy methods. The transition 
pressure provides a sensitive measure of small relative energy 
changes between the two phases (one a semiconductor and the 
other a semimetal). Experimentally, the transition pressure is 
well characterized. Density-functional theory (DFT) results 
have been unsatisfactory, exhibiting sensitivity to the particu- 
lar form of the exchange-correlation (xc) functional. Even the 
generally much more accurate diffusion Monte Carlo (DMC) 
metho d 2 ' 3 ' 4,5,6 has shown^ a noticeable fixed-node^ error. 

The phaseless auxiliary-field (AF) quantum Monte Carlo 
(QMC) AFQMC metho d 9 ' 10 : 11 provides a new alternative for 
ab initio many-body calculations to address electron corre- 
lation effects. All stochastic QMC method s 4 ' 9 ' 12 - 13 use pro- 
jection from a reference many-body wave function. In prin- 
ciple these methods are exact. In practice, however, the 
fermionic sign proble m 4 - 9 ' 14 - 15 ' 16 causes exponential growth of 
the variance with system size and projection time. Transient 
methods j 12 ' 17 ' 18 which maintain exactness while enduring the 
sign problem, can be very useful ;/ sufficiently accurate infor- 
mation can be obtained with a relatively short projection, as 
we illustrate in the present paper (Sec. lIV Al l. In general, how- 
ever, the sign problem must be completely eliminated (usually 
with an approximation) to achieve a general, efficient method 
for realistic systems. The majority of QMC calculations in 
fermion systems have been done in this form, for example 
with the fixed-node approximation 4 ^ in DMC, which has been 



the most commonly applied QMC method in electronic struc- 
ture. 

The phaseless AFQMC controls the sign problem with a 
global phase condition in the over-complete manifold of Slater 
determinants (in which antisymmetry is imposed). Since 
the antisymmetry ensures that each walker is automatically 
"fermionic", the tendency for the walker population to col- 
lapse to a global bosonic state is eliminated in this approach. 
It is reasonable to expect that an overall phase constraint 
applied in this manifold to be less restrictive^. Applica- 
tions indicate that this often is the case. In a variety of sys- 
tems AFQMC has demonstrated accuracy equaling or sur- 
passing the most accurate (non-exponential scaling) many- 
body computational methods. These include first- and second- 
row molecules^ *!?, transition metal oxide molecules, 2 ° simple 
solidsri^i post-d elements^ van der Waals systems^ molec- 
ular excited states, 24 and in molecules in which bonds are be- 
ing stretched or brokeni 24 ' 25 ' 26 Most of these calculations used 
a mean-field single determinant taken directly from DFT or 
Hartree-Fock (HE) for the trial wave function in the phase- 
less constraint. As a result, the phaseless AFQMC method 
reduces the reliance of QMC on the quality of the trial wave 
functioni 10 - 25 : 26 This is desirable in order to make QMC more 
of a general and "blackbox" approach. 

The use of a basis set is a second feature that distinguishes 
the AFQMC method from the standard DMC method ; 2 : 3 - 4 - 5 : 6 
The latter works in electron coordinate space. As a result, 
there is no finite basis set error per se in DMC. There are 
presently two main flavors of the phaseless AFQMC method, 
corresponding to two different choices of the one-electron 
basis: (i) planewave with norm-conserving pseudopotential 
(as widely adopted in solid state physics)) 9 ^ and (ii) Gaus- 
sian type basis sets (the standard in quantum chemistry).— In 
planewave AFQMC, convergence to the basis set limit is eas- 
ily controlled, as in DFT calculations, using the plane wave 
cutoff energy E cut . 

In this paper, planewave AFQMC is used to calculate the 
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relative energy differences between the two phases. The goal 
is to examine the accuracy of phaseless AFQMC, benchmark- 
ing the energy difference at the transition volumes against ex- 
periment and DMC results, and against exact free-projection 
AFQMC using smaller primitive cells. In the phaseless 
AFQMC approach, all but the error from the phaseless con- 
straint can be controlled systematically and driven essentially 
to zero. Comparison with exact AFQMC free-projection 
shows that the systematic error in the total energies due to the 
phaseless constraint is well within O-SnxEh/atom. Consistent 
with these internal benchmarks, the transition pressure calcu- 
lated from the phaseless AFQMC in large supercells is found 
to be in very good agreement with experiment. 

The paper is organized as follows. Several aspects of the 
AFQMC method, including the hybrid formulation and the 
reduction of weight fluctuation, are described in Sec. [II] This 
is followed by specific planewave AFQMC calculational de- 
tails in Sec. [Ill] Calculated results are presented and discussed 
in Sec. [IV] Finally, we summarize and conclude in Sec.[Vl 



Using Eq. (O effectively maps the two-body interaction onto a 
fictitious non-interacting Hamiltonian with coupling to auxil- 
iary classical fields {<Ji} = er. The operation of the one-body 
projector on a Slater determinant \ <fr) simply yields another de- 
terminant: \4>') = e v ^°" v |0). If |\|/ T ) in Eq. (Tj is expressed 
as a sum of Slater determinants (e.g., just one if |\I> T ) is a HF 
or DFT solution), the integral in Eq. (0]i can then be evaluated 
using Monte Carlo sampling over random walker streams.^! 

As discussed further in Sec. Ill Bl it is advantageous com- 
putationally to rewrite the two-body potential in Eq. ([5]), sub- 
tracting the mean-field contributio n 1 7 ' 1 8 i 28 i 32 prior to the HS 
transformation: 



V = — 



V m f) 2 + V • Vmf - 2 V mf 



(5) 



where v m f is generally chosen to be the expectation value of 
the v operator with respect to the trial wave function 



Vmf 
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(6) 



II. AFQMC METHODOLOGY 

This section reviews aspects of the AFQMC method in 
some detail. This is done to facilitate the discussion of sys- 
tematic errors in Sees. UII1 and UVl and to provide additional 
details on some phaseless AFQMC variants which are used 
in this paper. More complete descriptions of the phaseless 
AFQMC method can be found in Refs-ffiMDHisl 



A. AFQMC projection by random walks 

The ground state of a many -body system, is obtained 
by means of iterative projection from a trial wave function 
|* T ): 
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(1) 



where H = K + V is the Hamiltonian of the system, con- 
sisting of all one-body terms, K, and two-body terms, V. 
AFQMC implements the ground-state projection as random 
walks in the space of Slater determinants. The Trotter-Suzuki 
breakup 
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is used to separate the one- and two-body terms. Expressing 
V as a sum of the squares of one-body operators {vi} = v: 



V 



1- - 



(3) 



the Hubbard-Stratonovich (HS) transformation^^ is then 
used to express the two-body projector as a multidimensional 
integral 



n 



da, 



e -* 2 j2 e ^Fa % v % (4) 



B. Phaseless AFQMC 

In principle, the procedure in Eqs. (Q][4]) yields the ex- 
act ground state. The basic idea can be efficiently real- 
ized by branching random walks, as is used in Sec. IIV Al to 
carry out exact free-projection. In practice, however, a phase 
problem appears, because the repulsive Coulomb interaction 
gives rise to imaginary v, complex walkers \<f>), and complex 
(^ T \4>) overlaps, causing the variance to grow exponentially 
and swamp the signal. To control this problem, importance 
sampling and a phaseless approximation^ were introduced, 
yielding a stable stochastic simulation. The importance sam- 
pling transformation leads to a representation of the ground- 
state wave function as a weighted sum of Slater determinants 



l*o> = £ 



W0 
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A force bias term results in Eq. ©: 

da. 



ry =n 



'2-k 



(7) 
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The corresponding importance-sampled one-body propagator 
then takes the form 



da-g(cr)B{(T - &)W(a,&) 



(9) 



where {cr,} = <x, g(cr) is the multidimensional Gaussian 
probability density function with zero mean and unit width, 
and 



B(a -&) = e -rK/2 e Vf e -rK/2 ^ 

\<t>>) = B{* - a)\<j>) , 



(10) 

(11) 

(12) 
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The one-body operator B generates the random walker 
stream, transforming \<ft) into \<j)'), while W updates the 
weight factor w</,. 

The optimal choice of er, which cancels the weight fluctua- 
tion to O (y/r), is given by 



(13) 



where (ft is the determinant being propagated. Using this 
choice, the weight update factor W can be written as 9 



W(tr,a) 



e -r(* T |ff|0)/(* T |0) = e -rE L 



(14) 



where E L [<ft] is referred to as the "local energy" of \<ft). In 
practice, we use the average of two local energies to update 
the weight: 



W(tr,&) 



-t(B l [0']+£ l [0])/2 



(15) 



The total energy can be calculated using the mixed-estimate 
form, which is not variational 9 . 

The key to controlling the phase problem is to prevent a 
two-dimensional random walk in the complex (^ T \cft) -plane, 
thus avoiding the growth of a finite density at the origin. To 
do this, the phase rotation of the walker \<j>) is defined by 



A8 = arg 



(16) 



and the walker weight is "projected" to its real, positive value: 
jcos(A9)\W(cr,a-)\w 4> , \A6\ < ir/2 



0. 



otherwise 



(17) 



If the mean-field background is non-zero, its subtraction in 
Eq. (f5]l can lead to a reduction in the average rotation angle 
AO (and variance of the energy)i22i22 



C. AFQMC in hybrid form 

Most applications to date have used the phaseless AFQMC 
local energy formalism, described above. In planewave 
AFQMC, evaluating E L scales as O (N 2 M\ogM), while 
the propagation step [Eq. (O] scales as O (NM log M), using 
fast Fourier transforms.— Computation of the overlap matrix 
and other operations scale no worse than O (7V 2 M) . 

To reduce the frequency of evaluating E L , the most costly 
part of the calculation, we can use an alternative formulation, 
the "hybrid" form2i2! of the walker weight in Eq. (1121) . In 
the hybrid variant, only measurement evaluations of E L are 
needed. Since the autocorrelation time is typically 50-100 
times the time step r, this variant may be more efficient. The 
hybrid method tends to have larger variance than the local en- 
ergy method, however. The latter satistifies zero-variance in 
the limit of an exact \ ^>t), explicitly canceling out some O (r) 
terms. The two methods also have different Trotter behaviors, 
as illustrated in Sec. IIIIBI but they approach the same answer 
as r — > 0. The hybrid method is used for the large supercell 
calculations reported in this paper. 



D. Random walk bounds: controlling rare event fluctuations 

For any finite population of walkers, the stochastic nature of 
the simulation does not preclude rare events, which cause ex- 
tremely large population fluctuations. For example, a walker 
near the origin of the (4' T |(/))-plane can acquire a very large 
weight in a move \<ft') «— B\<ft) [Eq. (0], due to the occurence 
of a very large ($> T \(ft') / (^ T \(ft) ratio [Eq. (|T2j or Eq. (IT4b1. 
To circumvent the problem in a simulation of finite popula- 
tion, we apply a bound condition in the local energy method: 



(E° L -AE L )<E L ^}<(E° L + AE L ), 
where the width of the energy range AEl is defined as 

2 



AE L = \l- 



(18) 



(19) 



and where the average local energy value E^ is obtained by 
averaging E L measurements during the growth phased If E L 
goes outside this range, it is capped at the maximum or min- 
imum of the range. For a typical t (~ 0.05-E^ -1 ), the energy 
range allowed by Eq. ( fT8l is large (~ 12Eh), so E L is capped 
only in very rare instances. 

Similar bounds are introduced in the hybrid AFQMC 
method. Defining the hybrid energy as [compare Eqs. (TT2l 
and dm 



log W(a, <x) 



log 



<j ■ <j a ■ er 

2 



(20) 



the value of E H is bounded as 



(E° H - AE K ) < E H [cft} < (E° H + AE H ) , (21) 

where E^ is estimated as in Eq. ( fT8l . 

In addition, the walker weights are also bounded such that 
w<j> < Wmax at all times for a reasonable w max (typically set 
to the smaller of 100 or 0.1 times the size of the population). 
This bound is rarely triggered when the E s or E L bounding 
scheme is in place. 

Finally, a force-bias bound is applied in both the local en- 
ergy and hybrid methods. This prevents large modification 
of the orbitals when the denominator (^l^) in Eq. ( fT3l is 
small: 



1^1 < 1-0. 



(22) 



This bound is implicitly r-dependent, as seen in Eq. (1131) . We 
have found that the energy cap (E L or E u ) had the most effect 
in controlling weight fluctuations. 

It is important to note that the bounds being applied, while 
ad hoc, have well-defined limiting behavior. As r — > 0, the 
bounds on the physical quantities E^ and (v) both approach 
oo. The bounds only affect the Trotter error at finite r, but not 
the final answer when r is extrapolated to zero. 
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E. Exact calculations: unconstrained AFQMC 

To estimate the accuracy of phaseless AFQMC, calcula- 
tions using exact unconstrained "free" projection were carried 
out (Sec. lIHI) . In free projection, the weights { w$} are allowed 
to acquire a phase. This is implemented using a modified form 
of the hybrid method, where the mean-field average of the v 
operators is used as the force bias [instead of Eq. (1131) 1. 



>T V m f 



(23) 



This choice is equivalent to the subtraction of mean-field con- 
tribution to the two-body potential described in Eq. (0. The 
use of the mean-field background subtraction is essential in 
prolonging the stability of the simulation before the signal is 
lost to the phase problem. None of the bounds in the preceding 
subsection is applied in the free-projection calculations. 



III. AFQMC COMPUTATIONAL DETAILS FOR SILICON 
DIAMOND AND /3-TIN 

The present calculations are carried out with planewave 
based AFQMC (PW-AFQMC), which uses norm-conserving 
and separable Kleinman-Bylander 14 pseudopotentials to 
achieve efficient O (N 2 M log M) system size scaling, 
similar to planewave based DFT (PW-DFT) calculations. 
We first describe specific computational details of the 
planewave AFQMC calculations, including the pseudopoten- 
tial, planewave cutoff, and supercells. 

Convergence to the basis set limit is easily controlled, as in 
DFT calculations, using the plane wave cutoff energy E cut . 
Our calculations used E cut = 12.5i?h, which is the design 
cutoff of our Si pseudopotential (see below). For material sys- 
tems such as silicon, we have previously shown 11 that a good 
E cut at the DFT level, as determined by the norm-conserving 
pseudopotential, is sufficient to converge the two-particle cor- 
relations in AFQMC to within typical ststistical errors. In 
DFT calculations with the local density approximation (LDA), 
the total absolute energies of the diamond primitive cell using 
this E'cut has an error of ~ 0.43 mEh (as verified by using in- 
creasingly larger values of -Ecut)- Basis set convergence errors 
of energy differences are much smaller, of course. 

AFQMC calculations for large 54-atom 3x3x3 diamond 
and /3-tin supercells were done to obtain the transition pres- 
sures, after finite-size corrections, discussed below. Test cal- 
culations, such as pseudopotential tests and comparisons with 
benchmark exact AFQMC, were carried out for the smaller 
2-atom primitive unit cells. 

For each supercell and k-point, corresponding trial wave 
functions \^> T ) were taken as generated from DFT-LDA, us- 
ing the ABINIT code^. In /3-tin, random k points are used, 
rather than special points such as Monkhorst Pack sets, to re- 
move open-shell effects. For each k, our single-determinant 
trial wave function is thus unique and non-degenerate at the 
"Fermi surface." 

In the following subsections, aspects of the Si OPIUM pseu- 
dopotential are first discussed. The quality of the pseudopo- 
tential is assessed by comparing the equation of state (EOS) 



for the diamond and /3-tin structures with all-electron results 
within the framework of DFT. Next, efficient finite-size cor- 
rections are described, separately analyzing one-body errors, 
which are analogous to k sampling in PW-DFT, and two-body 
Coulomb finite size errors. 
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FIG. 1 : (Color online) Equation of state for the Si diamond and /3- 
tin phases, comparing all-electron (solid lines) and pseudopotential 
(dashed lines) DFT-LDA results. The inset shows the Gibbs energy 
as a function of the pressure. 



TABLE I: Quantities of the diamond and /3-tin phases of silicon 
computed with DFT, using all-electron LAPW and planewave pseu- 
dopotential methods. The xc functional used is the Perdew-Wang 
LDA. Volumes and lattice constants are expressed in atomic units 
(do and ao, respectively); energies are in eV; bulk moduli and pres- 
sures are in GPa. 



Quantity 


Pseudopotential 


LAPW 


diamond phase 






Equilibrium volume 


263.888 


266.474 


Equilibrium lattice constant 


10.182 


10.215 


Bulk modulus 


95.380 


95.327 


Cohesive energy 


5.413 


5.409 


/3-tin phase 






Equilibrium volume 


199.132 


200.364 


Bulk modulus 


114.760 


114.947 


Transition pressure 


7.67 


6.86 



A. Si pseudopotential quality 

The optimized design method^ was used to generate the 
Si pseudopotenital with OPIUM. — The atomic reference state 
was [Ne] 3s 2 3p 1,5 3d - 4 . All angular momentum channels 
(I = 0,1,2) used a cutoff radius r c = 2.08 Bohr, with 
/ = 2 as the local potential. The optimized design pseudo- 
wavefunction was expanded using five spherical Bessel func- 
tions with wave vector q c = 5.0 Bohr -1 , which corresponds 
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to a design E cut = 12.5 Eh, with a predicted plane wave con- 
vergence error of 1 mEh/atom for the absolute total energy. 
Explicit tests with DFT-LDA indicated errors several times 
smaller (see above), in both phases. 

To test the quality of the pseudopotential, the EOS for the 
diamond and /3-tin structures was compared to all-electron 
results within the framework of DFT. The results are shown 
in Fig. [TJ All-electron calculations were done using the 
ELk2£ full-potential LAPW program, and pseudopotential cal- 
culations with the planewave based ABINIT»- code (using the 
same OPIUM pseudopotential as in AFQMC). The DFT-LDA 
Perdew-Wang^S functional was used. In the all-electron and 
pseudopotential calculations, identical dense fc-point grids 
were used (6 x 6 x 6 in diamond and 16 x 16 x 16 in /3-tin). 
A temperature broadening of kgT = 0.05 eV was used in the 
/3-tin structure. Birch-Murnaghar42 fits were used to plot the 
Gibbs free energy. The agreement for the EOS between the 
pseudopotential and all-electron calculations is good, includ- 
ing the transition pressure values, which differ by ~ 0.8 GPa. 
These results are quantified in Table [I] 



B. Trotter errors 

The transition pressure calculations in Sec. IIVI were done 
for 3x3x3 supercells, using a Trotter time step of 
t = 0.025 E^ 1 . In benchmarking exact AFQMC results 
in Sec. |IV Al extraplotation to r — > was examined care- 
fully for lxlxl primitive cells for the phaseless local- 
energy and hybrid AFQMC methods as well as for exact free- 
projection. Not surprisingly, extrapolation errors largely can- 
cel between the two structures. For example, the residual er- 
rors at t = 0.025-E^ 1 are 1.7(1) and 1.5(1) mE h for diamond 
and /3-tin primitive cells, respectively. 

We also did several tests at larger supercell sizes. The resid- 



ual error at r = 0.025£' h " of a 3 x 3 x 3 diamond structure 
supercell was estimated to be 1.6(4) mEh (normalized to the 
primitive cell), very similar to the value of 1.7(1) mEh for the 
corresponding lxlxl primitive cell. No explicit Trotter cor- 
rections were applied, therefore, in calculating the transition 
pressure, given the error cancellation between the two struc- 
tures and the fact that the estimated residual errors in even 
the absolute energy are not significantly larger than the QMC 
statistical errors. 



and converges slowly, largely because two-body interactions 
are long-ranged, causing FS effects to persist to large system 
sizes. Alternatively, FS correction schemes can be usedi 41 ' 42 ' 43 
Both one- and two-body FS corrections 4 '' 43 must be applied 
to achieve efficient convergence. One-body effects are related 
to BZ k-point sampling. These can be largely corrected, using 
DFT calculations to estimate quadrature errors. In metals such 
at /3-tin, BZ intergration errors are aggravated by open-shell 
effects. Twist-averaged boundary conditions 43,44 can be used 
in this case to further reduce residual one-body errors, as is 
done here for the /3-tin phase. The one-body FS correction is 
given byi 3 - 



AG 



shell 



E L 



Et 



(24) 



namely by subtracting the DFT energy at the same k vector 
(^dfT) an( j a( j(jing the DFT energy obtained with a dense 
k grid (E DFT ). Figure [2] shows the reduced variation of 
the AFQMC total energy after this correction is applied, for 
3x3x3 /3-tin supercell. Averaging over the 9 randomly cho- 
sen k points before the correction results in a statistical error 
(combined error of the nine random data points each of which 
has a statistical error bar) of 1 mEh, while averaging after the 
correction reduces the combined error to 0.6 mEh. As men- 
tioned, random k points rather than special points were used 
to remove open-shell effects in metals and ensure that the trial 
wave function is non-degenerate. 




C. Finite-size errors 

Independent-particle methods, such as DFT or HF, can use 
Bloch's theorem to perform calculations in crystals, using 
only the primitive unit cell. The macroscopic limit is achieved 
by k-point quadrature in the Brillouin zone (BZ). Many-body 
methods, by contrast, must be performed for individual su- 
percells. The resulting finite-size (FS) errors often can be 
more significant than statistical and other systematic errors. 
Eliminating or reducing the FS errors is crucial, therefore, 
to achieve accurate results. The brute force extrapolations 
approach, using increasingly larger supercells, is expensive 



FIG. 2: (Color online) AFQMC and LDA Si /3-tin energies at nine 
randomly-chosen k points in the Brillouin zone of 3 x 3 x 3 supercell. 
The AFQMC and LDA one-body FS errors are correlated, and the 
correction in Eq. J24t reduces the variation in the QMC data. The 
two-body FS error is significant even with a 54-atom supercell, but 
is essentially independent of k-points. 



The two-body FS error comes from the artificially induced 
periodicity of the long-range electron-electron Coulomb re- 
pulsion, due to the use of periodic boundary conditions. This 
error can be reduced significantly, using the post-processing 
correction scheme of Kwee et al. which is based on a finite- 
size DFT xc functional, corresponding to the finite-sized su- 
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percell. This two-body correction is given by 

A £Coulomb = ^DFT.oo _ ^DFT.L (25) 

where E^' 00 [= E^ 7 in Eq. (23)] is the DFT energy 
computed with the usual LDA xc functional, while E^ FT L 
is the DFT energy computed with the finite-size LDA xc 
functional^ The k-dependence of /\E^° alomh is very small 
compared to that of the one-body correction shown in Fig. 12 
with variations of ~ O.lm-Eh in /3-tin. 

The total FS correction is the result of applying the one- 
and two-body correction terms, Eqs. d24i i and d25l ). respec- 
tively. This is of course equivalent to applying AE^ = 
e dft _ ^DFT-i to ^ raw AFQMC energies. The corrected 
energies are averaged over the k points. The net effect of ap- 
plying both FS corrections is to decrease the energy difference 
at the transition volumes from 34(1) to 29(1) mE^. With these 
combined FS corrections, the residual errors in the absolute 
energies from 3x3x3 supercells are expected to be small in 
silicon.- 4 - Error cancellation in the energy difference between 
/3-tin and diamond structures further reduces the error in the 
calculated transition pressure. 

IV. RESULTS AND DISCUSSION 

A. Benchmarking the phaseless approximation with exact 
free-projection AFQMC 

The fermionic sign/phase constraints used by QMC meth- 
ods generally introduce uncontrolled approximations. Exam- 
ples include the DMC fixed-node approximation and AFQMC 
phaseless constraint. Except where benchmarks with exact 
methods or experiment are available for comparison, the cor- 
responding constraint errors are difficult to quantify. In this 
section, we show that exact free-projection calculations are 
feasible for the primitive diamond and /3-tin structures, using 
planewave AFQMC on a large parallel computing platform. 
Comparison with the corresponding approximate phaseless 
AFQMC calculations shows that the systematic error due to 
the phaseless constraint is small (within 0.5 m£h/atom), as 
described below. 

As illustrated in Fig. [3] for the diamond structure, free- 
projection to the ground state can be achieved in the primi- 
tive cell using large walker populations. The free-projection 
calculation was done with a target population size of two mil- 
lion walkers, using about 2000 cores at the NCCS Jaguar 
XT4 computer at Oak Ridge National Laboratory. An ac- 
ceptable signal-to-noise ratio is sustained for sufficiently long 
imaginary times. For projection times /3 >~ 16E^~ , how- 
ever, growing fluctuations, due to the phase problem, begin to 
emerge. Eventually the fluctuations become severe enough to 
destroy the Monte Carlo signal. 9 The energy measurement for 
this benchmark is taken after the walkers are sufficiently equi- 
librated, (3 > 10E^ . Similar calculations were performed in 
the /3-tin structure. 

Figures |4j and [5] display extrapolations of the Trotter time 
step, r, for the diamond and /3-tin structures, respectively. The 
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FIG. 3: (Color online) Comparison of phaseless AFQMC with ex- 
act free-projection. The calculations shown are for the primitive cell 
in the diamond structure at the experimental equilibrium lattice con- 
stant (10.264 ao). The calculations used L as the k-point, and a time 
step size of r = 0.025-Ef 1 . 

figures compare the local-energy and hybrid phaseless meth- 
ods with free-projection results. As r — > 0, the local-energy 
and hybrid phaseless methods are seen to converge to the same 
result, as expected, but with different slopes. To leading or- 
der, free-projection shows t 2 behavior, while the local-energy 
and hybrid methods have linear r behavior, since the phase- 
less constraint (and the bounds in Sec. Ill Db in the latter two 
methods break the quadratic scaling in Eq. (0. The hybrid 
method in Figs.|4]and|5]is seen to have the largest slope. The 
Trotter behaviors of the respective methods are similar in the 
diamond and /3-tin structures. 

The error in the total energy caused by the phaseless ap- 
proximation, after extrapolation to r — > 0, is about 0.7mi?h 
(or 0.35mEh per atom) for diamond and 0.8m£"h in /3-tin. 
Note that the energy calculated from the phaseless approxi- 
mation using the mixed-estimate [Eq. ( fl~4b l is not variational 9 . 
Indeed in both cases above it is below the exact result. 

TABLE II: AFQMC energies of Si diamond and /3-tin phases for 
the 2-atom primitive cells, with volumes of 40.07 A 3 and 30.00 A 3 , 
respectively. The /3-tin c/a ratio was set to 0.552. Phaseless AFQMC 
values are from the local energy formalism. Energies are in Eh units. 



phase 


reduced k vector 


free projection 


phaseless 


diamond 


(0.5,0.5,0.5) 


-8.05485(8) 


-8.0555(1) 


/3-tin 


(0.25,0.25,0.25) 


-8.1256(1) 


-8.1264(1) 


/3-tin 


(0.25,-0.25,0.25) 


-8.0023(2) 


-8.0017(2) 



In Tablelnl we list the absolute energies for three cases from 
our free-projection calculations, which can serve as bench- 
marks in the future. All the energies have been extrapolated 
to r — > 0. The corresponding phaseless results are also listed, 
which are in very good agreement with the exact results. We 
note that the two twist boundary conditions (k points) in /3- 
tin show different behaviors. In k = (0.25, 0.25, 0.25) the 
phaseless energy is below the exact value (as in the diamond 
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FIG. 4: (Color online) Trotter time step r extrapolation for the di- 
amond structure, as in Fig. [3] comparing local-energy and hybrid 
phaseless methods with free-projection. 
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/3-tin structures. For /3-tin, we used the experimental value of 
c/a = 0.5504& (It was shown in Ref. |7|that the dependence 
on c/a is weak.) Twist averaged boundary conditions used the 
single special k-point of Baldereschiii for the diamond struc- 
ture, and nine random k points for the /3-tin phase. The total 
energies lead to a "raw" transition pressure of 15.1(3) GPa. 

To compare with experiment, corrections are required to ac- 
count for zero-point motion and thermal effects. We apply 
these corrections as given in Ref. 0: 1) a zero-point motion 
lowering of 0.3 GPa; 2) a room-temperature quasi-harmonic 
estimate of the relative stabilization of the /3-tin phase, which 
lowers the pressure by 1.15 GPa. This would give a tran- 
sition pressure of 13.6(3) GPa. In addition, standard mean- 
field pseudopotentials generated from LDA or HF, such as 
the one used in the present paper, do not account for man 
body effects in the core. A correction was estimated in Ref. 
by explicitly including a many-body core polarization poten- 
tial (CPP), which further lowers the pressure by ~ 1.2 GPa. 
Assuming that our LDA pseudopotential is similar to that in 
Ref. we apply the same correction. Table [III] reports our 
final result, and compares it to experiment and to other theo- 
retical results. Corrections (1) and (2) have also been applied 
to the DFT results. 



TABLE III: The transition pressure P t of the diamond to /3-tin tran- 
sition in silicon. The AFQMC result is listed together with exper- 
imental and other theoretical results. To compare with experiment, 
appropriate corrections have been applied to the theoretical results 
(see text). 



Method 



Pt (GPa) 



Imaginary time step X (E ) 



LDA^ 
GGA (BP)^ 
GGA (PW91)^ 
DMGi 
AFQMC 
Experiment 



6.7 

13.3 

10.9 

16.5(5) 

12.6(3) 

10.3-12.5 



FIG. 5: (Color online) Trotter time step r extrapolation for the /3-tin 
structure, comparing the local-energy phaseless method with free- 
projection. Calculations are for the primitive unit cell with volume 
of 30 A 3 and c/a = 0.552 at the reduced k = (0.25, 0.25, 0.25). 



case), while in the other, the phaseless energy is above. This 
manifests the varying quality of the trial wave function at the 
different k-points, which shift the single-particle energy lev- 
els differently. Also, the FS effects are clearly very large in 
these small cells, where the energy for the first k-point is in 
fact below that of the diamond result. 



B. Transition pressure 

AFQMC calculations were done for 3x3x3 supercells 
(containing 54 silicon atoms), at the experimental transition 
volumes,^ 36.30 A 3 for the diamond and 27.91 A 3 for the 



The transition pressure is not very sensitive to the choice 
of transition volumes. For example, using the DMC predicted 
volumes instead of the experimental values changes the en- 
ergy difference by only 0.01 eV, from 0.49 to 0.50 eV, reduc- 
ing the transition pressure by less than 0.3 GPa. 

The best calculation to date with the highest level of theory 
is the DMC calculations in Ref. 7. Compared to experiment, 
the somewhat overestimated DMC Pt = 16.5(5) value was 
attributed to the fixed-node error^ This seems consistent with 
our results. The DMC discrepancy corresponds to a larger 
"raw" energy difference of 19.2 mE^ between the two phases, 
compared to 15.1(3) mEh for phaseless AFQMC. As shown 
in the previous subsection, the error due to the phaseless ap- 
proximation (~ 0.5 mEh) appears to be an order of magnitude 
smaller than this. Our calculations show that experiment and 
theory are in quantitative agreement on the diamond to /3-tin 
transition. 
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V. SUMMARY 

We have applid the phaseless auxiliary-field quantum 
Monte Carlo method to study the pressure-induced structural 
phase transition from diamond to /3-tin in silicon. This is 
a recently developed non-perturbative, many-body approach 
which recovers electron correlation by explicitly summing 
over fluctuating mean-field solutions with Monte Carlo. The 
only source of error which can not be systematically driven to 
zero is that of the global phase constraint, used to control the 
sign/phase problem. We quantified the systematic error from 
this phaseless approximation by exact unconstrained AFQMC 
calculations in the primitive cell, carried out on large parallel 
computers. In both structural phases the error was found to 
be well within 0.5mEh/atom. A transition pressure was cal- 
culated form the energy difference between the two phases 
at the experimental transition pressure, using 54-atom super- 
cells. Twist-averaging boundary condition and finite-size cor- 
rections were applied, which greatly accelerates the conver- 
gence to the thermodynamic limit. After corrections for zero- 
point effect, thermal effect, and the (lack of) core-polarization 
in the pseudopotential, the AFQMC results yield a transition 
pressure of 12.6 ± 0.3 GPa, compared to experimental values 



of 10.3-12. 5GPa. 

The good agreement between the phaseless AFQMC result 
and experiment is consistent with the internal benchmark with 
unconstrained AFQMC. Our analysis indicates that the pos- 
sible combined error from the calculations should be below 
1 GPa. These include pseudopotential transferability errors 
and core-polarization effect, residual finite-size errors, and the 
error from the phaseless approximation. 
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